% Figures for the paper
% 2017-2-27
% Tentative plan: We provide figures for houston
clc; clear all; close all;
workpath00 = pwd;

addpath('toolbox_xt');
addpath('toolbox_plot');

% predX_M1: covariates at the tract centroid is also estimated via
% hierachical model 
% M1: includes intercept and log(1+dist), coefficients are varying with Area
% M2: includes intercept and log(1+dist) and log(dist 2 coast), coefficients are varying with Area
estXmode = 3; %1=M1, 2=M2, 3=M3, estXmode 3 is the one we used in the paper

% !!!!
% ONLY WORKS FOR TABLE (NOT FIGURES)
% NOW FIGURES ARE UPDATED TOO (5/25/2017)

%% Figure option
cityid = 520; %520 for Atlanta, 1600 for Chicago, 4480 LA, 5600 NY, Houston 3360
nyear = 6;
modelind = 6; %2: linear, 6: quartic, size-cubic
islog  = 2; %2 for log(1+x)

do_tab   = 1; %land value table

%% Load Data and preliminary stuff
load('data_all');
% load('results_full_logp1.mat');
% load('results_full_logp1_A.mat');
load('results_full_logp1_A_s3_REML.mat');
load('results_estX_REML.mat');
if estXmode == 1
    predX = estX_M1.predX;
    predM = estX_M1;
elseif estXmode == 2
    predX = estX_M2.predX;
    predM = estX_M2;
elseif estXmode == 3
    predX = estX_M3.predX;
    predM = estX_M3;
end

% other info
ncity = size(ZZ2,1);

% Spatial function
ngrid = 100;
xgrid = linspace(0,60,ngrid)';
logxgrid = log(1+xgrid);

logxgrid1 = log(1+xgrid);
logxgrid4 = [logxgrid1, logxgrid1.^2, logxgrid1.^3, logxgrid1.^4];



% for i=2:1:5
%     eval(['temp_m = est_s',num2str(i),';']);
%     eval(['sp',num2str(i),'= get_sp2(temp_m,islog);']);
%     eval(['sp',num2str(i),'.sp     = sp',num2str(i),'.spfun(xgrid);']);
%     eval(['sp',num2str(i),'.xgrid  = xgrid;']);
%     eval(['est_s',num2str(i), '.sp = sp',num2str(i),';']);
% end



%% Table: Computation of the land value
if do_tab == 1
    temp_x = XX;
    temp_y = YY;
    
    eval(['temp_est = est_s',num2str(modelind),';']);
    %m_XX = mean(XX);
    %m_XX2 = mean([XX(:,14).^2, XX(:,14).^3]); %nonlinear term in size
    
    % CITY SPECIFIC COVARIATES
    [m, big_m_XX, info] = xt_mean(XX,xtflag,xtval);
    [m, big_m_XX2, info] = xt_mean([XX(:,14).^2, XX(:,14).^3],xtflag,xtval);
    XX2 = [XX(:,14).^2, XX(:,14).^3];
    
    npol = size(temp_est.post.m,1)-6; %number of polynomials
    temp_min_dist= log(1+ tract.min_dist);
    grid = [temp_min_dist, temp_min_dist.^2, temp_min_dist.^3, temp_min_dist.^4];
    grid = grid(:,1:npol);
    
    
    tab_area = [];
    tab_ua   = [];
    tab_land = [];
    
    tab_lv = [];
    for j=1:1:ncity;
        
        % city specific predictors
        %m_XX = big_m_XX(j,:);
        %m_XX2 = big_m_XX2(j,:);
        
        temp_XX  = XX(xtflag == xtval(j),:); %covariates in this city
        temp_XX2 = XX2(xtflag == xtval(j),:);
        temp_DD  = DD(xtflag == xtval(j),:);
        temp_ZZ  = ZZ(xtflag == xtval(j),:);
        
        % temporary variables
        cityid = xtval(j);
        temp_sel  = tract.combofips==cityid;
        temp_grid = grid(temp_sel,:);
        temp_lon = tract.lon(temp_sel);
        temp_lat = tract.lat(temp_sel);
        temp_d2c = tract.d2c(temp_sel);
        
        temp_area = tract.landsqmi_ua(temp_sel); %urban area
        temp_area_n = temp_area/sum(temp_area); % normalized area
        
        temp_area2 = tract.landsqmi(temp_sel); %urban area
        temp_area2_n = temp_area2/sum(temp_area2); % normalized area
        
        temp_n = size(temp_grid,1);
        
        % contructing distance dependent covariates for prediction
        % predX: we estimate and predict outside of this m-file.
        temp_predX1 = predX(temp_sel,1:14);
        temp_predX2 = predX(temp_sel,15:16);
        
        
        % land value at tract center in time t
        % time-invariant part of the value
        temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
            + temp_est.c10*log(temp_d2c) + temp_est.post.sig2/2; % value does not depend on time
        
%         % log-normal correction for the predictor uncertainty
%         temp_est.betbet = [temp_est.bet; temp_est.bet2];
%         temp_val_nv = m_XXX*temp_est.betbet + temp_est.betbet'*cov_XXX*temp_est.betbet ...
%             + temp_est.c10*log(temp_d2c) + temp_est.post.sig2/2; % value does not depend on time
        
        temp_val = zeros(temp_n, nyear);
        for tind = 1:1:nyear
            tdum = zeros(1,nyear);
            tdum(tind) = 1; %this is for alp_jt
            tdum(1) = 1; %this is alp_j
            temp_C = [repmat(tdum, temp_n,1), temp_grid];
            
            temp_val(:,tind) = temp_val_nv + temp_C*temp_est.post.m(:,j) + 1/2*diag(temp_C*temp_est.post.V(:,:,j)*temp_C');
        end
        temp_val = exp(temp_val);
        
        % [1] Integrated using urban area, average over time
        temp_val_ua = sum(temp_val.*repmat(temp_area_n,1,nyear),1);
        temp_val1 = mean(temp_val_ua); %average over years
        
        % [2] Integrated using land area, average over time
        temp_val_land = sum(temp_val.*repmat(temp_area2_n,1,nyear),1);
        temp_val2 = mean(temp_val_land);
        
        % [3] Total using urban area, average over time
        temp_val_ua_tot = sum(temp_val.*repmat(temp_area*640,1,nyear),1); %acre correction
        temp_val3 = mean(temp_val_ua_tot); %average over years
        
        % [4] Total using land area, average over time
        temp_val_land_tot = sum(temp_val.*repmat(temp_area2*640,1,nyear),1); % acre correction
        temp_val4 = mean(temp_val_land_tot); %average over years
        
        
        % [5] 1, 5, 10, 20, 40 miles from the city hall
        temp_grid5 = log(1+[0,1/2,1,5,10,20,40])';
        temp_grid5 = [temp_grid5, temp_grid5.^2, temp_grid5.^3, temp_grid5.^4];
        temp_grid5 = temp_grid5(:,1:npol);
        
        mean_temp_d2c = repmat(mean(log(temp_d2c)), size(temp_grid5,1),1);
        
        % contructing distance dependent covariates for prediction
        n_temp_pred = size(temp_grid5,1);
        temp_pred_grid5 = temp_grid5(:,1);
        temp_pred_ZZ = repmat(temp_ZZ(1,2), n_temp_pred,1);
        
        if estXmode == 1
            new_x = [ones(n_temp_pred,1), temp_pred_ZZ, ...
                temp_pred_grid5, temp_pred_ZZ.*temp_pred_grid5];
            new_z = {[ones(n_temp_pred,1), temp_pred_grid5]};
        elseif estXmode ==2
            new_x = [ones(n_temp_pred,1), temp_pred_ZZ, ...
                temp_pred_grid5, temp_pred_ZZ.*temp_pred_grid5, ...
                mean_temp_d2c, temp_pred_ZZ.*mean_temp_d2c, ...
                ];
            new_z = {[ones(n_temp_pred,1), temp_pred_grid5, mean_temp_d2c]};
        elseif estXmode ==3
            new_x = [ones(n_temp_pred,1), temp_pred_ZZ, ...
                temp_pred_grid5, temp_pred_ZZ.*temp_pred_grid5, ...
                mean_temp_d2c, ...
                ];
            new_z = {[ones(n_temp_pred,1), temp_pred_grid5]};
        end
        new_g = repmat(cityid, n_temp_pred,1);
        
        pred_y = [];
        for  vind = 1:1:16;
            eval(['temp_m = predM.model', num2str(vind), ';']);
            new_y = predict(temp_m,new_x,new_z,new_g);
            pred_y = [pred_y, new_y];
        end
        temp_predX1 = pred_y(:,1:14);
        temp_predX2 = pred_y(:,15:16);
 
        
        % time-invariant part of the value
        temp_val_nv = temp_predX1*temp_est.bet + temp_predX2*temp_est.bet2 ...
            + temp_est.c10*mean_temp_d2c + temp_est.post.sig2/2; % value does not depend on time
        
        %temp_est.betbet = [temp_est.bet; temp_est.bet2];
        %temp_val_nv = m_XXX*temp_est.betbet + temp_est.betbet'*cov_XXX*temp_est.betbet ...
        %    + temp_est.c10*mean_temp_d2c + temp_est.post.sig2/2; % value does not depend on time
        
        temp_val5 = zeros(size(temp_grid5,1), nyear);
        for tind = 1:1:nyear
            tdum = zeros(1,nyear);
            tdum(tind) = 1; %this is for alp_jt
            tdum(1) = 1; %this is alp_j
            temp_C = [repmat(tdum, size(temp_grid5,1),1), temp_grid5];
            
            temp_val5(:,tind) = temp_val_nv + temp_C*temp_est.post.m(:,j) + 1/2*diag(temp_C*temp_est.post.V(:,:,j)*temp_C');
        end
        temp_val5 = mean(exp(temp_val5),2)'; %year average
        
        % [6] Integrated urban land value over time
        temp_val6 = temp_val_ua;
        
        % [0] Simple average
        temp_val0 = mean(exp(temp_y(xtflag == xtval(j),:)));
        
        % [00] Simple average (weigted)
        temp_a = exp(temp_y(xtflag == xtval(j),:));
        temp_b = exp(temp_x(xtflag == xtval(j),14)); %log(size)
        temp_val00 = sum( (temp_a.*temp_b)/sum(temp_b)); %weighted average
        
        % store
        temp_tab = [temp_val0,temp_val00, temp_val1, temp_val2, temp_val3, temp_val4, temp_val5, temp_val6];
        tab_lv = [tab_lv; temp_tab];
        
        tab_area = [tab_area; [sum(temp_area), sum(temp_area2)]]; %ua and land area
        tab_ua = [tab_ua; temp_val_ua];
        tab_land = [tab_land; temp_val_land];
        
    end
    tab_head = {'combofips', 'cmsa', 'pmsa', 'n', 'simple', 'simple (weighted)', 'LV (UA)', 'LV (Land)', 'LV-total (UA)', 'LV-total (Land)', ...
        'LV (center)', 'LV (1/2 mile)', 'LV (1 mile)', 'LV (5 mile)', 'LV (10 mile)', 'LV (20 mile)', 'LV (40 mile)', ...
        'LV-2005', 'LV-2006', 'LV-2007', 'LV-2008', 'LV-2009', 'LV-2010', ...
        };
    
    % combofips information
    tab_name = {};
    for j=1:1:ncity
        temp_sel = name.pmsafips==xtval(j);
        
        tab_name = [tab_name; [name.pmsafips(temp_sel,:), name.cmsaname(temp_sel,:), name.pmsaname(temp_sel,:), name.ntot(temp_sel,:)]];
        
    end
    
    big_tab = [tab_head; [tab_name, num2cell(tab_lv)]];
    
    % save
    tabpath = [workpath00, filesep, 'figure_draft_rev'];
    chk_dir(tabpath);
    
    savefilename = ['tab_lv_log',num2str(islog),'m',num2str(modelind), '_estX_M', num2str(estXmode)];
    cd(tabpath);
    xlswrite(savefilename, big_tab, 'land_values (MSAs)', 'A1');
    cd(workpath00);
    
    %% US average
    % [1] US-average per acre, weighted by urban area
    lv_us1 = sum(tab_ua.*repmat(tab_area(:,1)/sum(tab_area(:,1)),1,nyear),1); %average, per acre
    % [2] US-average per acre, weighted by land area
    lv_us2 = sum(tab_land.*repmat(tab_area(:,2)/sum(tab_area(:,2)),1,nyear),1); %average, per acre
    % [3] US-total ua
    lv_us3 = sum(tab_ua.*repmat(tab_area(:,1)*640,1,nyear),1); %acre correction
    % [4] US-total land
    lv_us4 = sum(tab_land.*repmat(tab_area(:,2)*640,1,nyear),1); %acre correction
    
    
    tab_us = [lv_us1', lv_us2', lv_us3', lv_us4'];
    tab_name = {'2005', '2006', '2007', '2008', '2009', '2010'}';
    tab_head = {'year', 'US-per-acre (ua)', 'US-per-acre (land)', 'US-total (ua)', 'US-total (land)'};
    big_tab2 = [tab_head; [tab_name, num2cell(tab_us)]];
    
    % save
    tabpath = [workpath00, filesep, 'figure_draft_rev'];
    chk_dir(tabpath);
    
    savefilename = ['tab_lv_log',num2str(islog),'m',num2str(modelind), '_estX_M', num2str(estXmode)];
    cd(tabpath);
    xlswrite(savefilename, big_tab2, 'land_values (US)', 'A1');
    cd(workpath00);
end




